/* ============================================================================
 * I B E X - Interval definition
 * ============================================================================
 * Copyright   : IMT Atlantique (FRANCE)
 * License     : This program can be distributed under the terms of the GNU LGPL.
 *               See the file COPYING.LESSER.
 *
 * Author(s)   : Gilles Chabert
 * Bug fixes   : Gilles Trombettoni, Bertrand Neveu
 * Created     : Dec 05, 2011
 * Last update : Oct 30, 2019
 * ---------------------------------------------------------------------------- */

#ifndef _IBEX_INTERVAL_H_
#define _IBEX_INTERVAL_H_

#include <array>
#include <math.h>
#include "ibex_Exception.h"

/* ========================================================*/
/* The following header file is automatically generated by
 * the compilation.
 */
#include "ibex_Setting.h"

#include "ibex_IntervalLibWrapper.h"

/** \brief NEG_INFINITY: <double> representation of -oo */
#define NEG_INFINITY IBEX_INTERVAL_LIB_NEG_INFINITY
/** \brief POS_INFINITY: <double> representation of +oo */
#define POS_INFINITY IBEX_INTERVAL_LIB_POS_INFINITY

#ifdef _MSC_VER
#include <intrin.h>
#define __builtin_popcount __popcnt
#define __builtin_powi(__x,__n) ((double)::pow((double)(__x),(int)(__n)))
#endif // _MSC_VER


namespace ibex {

class IntervalVector;
class IntervalMatrix;
class ExprConstant;

/** \defgroup arithmetic Interval Arithmetic */

/** \ingroup arithmetic */
/*@{*/

/**
 * \brief Sets the rounding direction mode of the FPU towards -oo.
 */
void fpu_round_down();

/**
 * \brief Sets the rounding direction mode of the FPU towards +oo.
 */
void fpu_round_up();

/**
 * \brief Sets the rounding direction mode of the nearest
 */
void fpu_round_near();

/*
 * Sets the rounding direction mode of the FPU towards zero.
 */
void fpu_round_zero();

/**
 * \brief Return the previous float
 */
double previous_float(double x);

/**
 * \brief Return the next float
 */
double next_float(double x);


/*@}*/

/**
 * \ingroup arithmetic
 *
 * \brief Interval
 *
 * This class defines the interval interface of IBEX and encapsulates an interval "itv" whose
 * type depends on the chosen implementation (currently: Gaol, Bias of filib).
 *
 * Note that some functions of the Gaol interval interface do not appear here (like "possibly relations")
 * because there are not used by ibex; while other have been introduced (like "ratio_delta"). Some
 * functions are also renamed to match more conventional use.
 *
 * Note that with filib several precision and mode are available. We choose :
 * base type = double
 * rounding_strategy = native_switched
 * interval_mode = i_mode_extended_flag
 *
 */
class Interval {
  public:
	/** \brief Create (-oo,+oo). */
    Interval();

    /** \brief Create [a,b].  */
    Interval(double a, double b);

    /** \brief Create [a,a]. */
    Interval(double a);

    /** \brief Create [a, a] from a fixed-sized array of size 1. */
    Interval(std::array<double, 1> array);

    /** \brief Create [a, b] from a fixed-sized array of size 2. */
    Interval(std::array<double, 2> array);

    /** \brief True iff *this and x are exactly the same intervals. */
    bool operator==(const Interval& x) const;

    /** \brief True iff *this and x are not exactly the same intervals. */
    bool operator!=(const Interval& x) const;

    /** \brief Set this interval to the empty set. */
    void set_empty();

    /** \brief Set *this to x.
     */
    Interval& operator=(const Interval& x);

    /** \brief Set *this to d.
     */
    Interval& operator=(double x);

    /** \brief Intersection of *this and x.
     * \param x - the interval to compute the intersection with.*/
    Interval& operator&=(const Interval& x);

    /** \brief Union of *this and x.
     * \param x - the interval to compute the hull with.*/
    Interval& operator|=(const Interval& x);

    /**
     * \brief Add [-rad,+rad] to *this.
     *
     * Return a reference to *this.
     */
    Interval& inflate(double rad);

	/**
	 * \brief Absolute and relative inflation.
	 *
	 * [x] <- mid[x] + delta*([x]-mid[x]) + chi*[-1,+1]
	 *
	 * \return *this.
	 */
	Interval& inflate(double delta, double chi);

    /** \brief Lower bound.
     *
     * Return the lower bound of *this. */
    double lb() const;

    /** \brief Upper bound.
     *
     * Return the upper bound of *this. */
    double ub() const;

    /** \brief Midpoint.
     *
     * Returns the midpoint of *this.
     * The return point is guaranteed to be included in *this
     * but not necessarily to be the closest floating point
     * from the real midpoint.
     *
     * Cases are:
     * - \emptyset  -> Quiet NaN
     * - [-oo, +oo] -> midP = 0.0
     * - [-oo, b]   -> midP = -MAXREAL
     * - [a, +oo]   -> midP = MAXREAL
     * - [a, b]     -> midP ~ a + .5*(b-a) */
    double mid() const;

    /**
     * \brief Radius.
     *
     * Return the diameter of *this.
     * By convention, 0 if *this is empty.*/
    double rad() const;

    /**
     * \brief Diameter.
     *
     * Return the diameter of *this.
     * By convention, 0 if *this is empty.*/
    double diam() const;

    /**
     * \brief Mignitude.
	 *
     * Returns the mignitude of *this:
     * <lu>
     * <li> +(lower bound)  if *this > 0
     * <li> -(upper bound) if *this < 0
     * <li> 0 otherwise.
     * </lu> */
    double mig() const;

    /**
     * \brief Magnitude.
	 *
     * Returns the magnitude of *this:
     * mag(*this)=max(|lower bound|, |upper bound|). */
    double mag() const;

    /**
     * \brief True iff this interval is a subset of \a x.
     *
     * \note Always return true if *this is empty.
     */
    bool is_subset(const Interval& x) const;

    /**
     * \brief True iff this interval is a subset of \a x and not \a x itself.
     *
     * \note In particular, (-oo,oo) is not a strict subset of (-oo,oo)
     * and the empty set is not a strict subset of the empty set although
     * in both cases, the first is inside the interior of the second.
     */
    bool is_strict_subset(const Interval& x) const;

    /**
     * \brief True iff this interval is in the interior of \a x.
     *
    * \note In particular, (-oo,oo) is in the interior of (-oo,oo)
     * and the empty set is in the interior of the empty set.
     * \note Always return true if *this is empty.
     */
    bool is_interior_subset(const Interval& x) const;

    /**
     * \brief True iff this interval is in the relative interior of \a x.
     *
     * When x is degenerated, the relative interior of x is x itself
     * (it is not an open set). Otherwise, the relative interior of x
     * is the interior (in the usual meaning).
     */
    bool is_relative_interior_subset(const Interval& x) const;

    /**
     * \brief True iff this interval is in the interior of \a x and different from x.
     *
     * \note In particular, (-oo,oo) is not "strictly" in the interior of (-oo,oo)
     * and the empty set is not "strictly" in the interior of the empty set.
     */
    bool is_strict_interior_subset(const Interval& x) const;

    /**
     * \brief True iff this interval is a superset of \a x.
     *
     * \note Always return true if x is empty.
     */
    bool is_superset(const Interval& x) const;

    /**
     * \brief True iff this interval is a superset of \a x different from x.
     *
     * \see #is_strict_subset(const Interval&) const.
     */
    bool is_strict_superset(const Interval& x) const;

    /**
     * \brief True iff *this contains \a d.
     *
     * \note d can also be an "open bound", i.e., infinity.
     * So this function is not restricted to a set-membership
     * interpretation.
     */
    bool contains(const double& d) const;

    /**
     * \brief True iff the interior of *this contains \a d.
     *
     */
    bool interior_contains(const double& d) const;

    /**
     * \brief True iff *this and \a x intersect.
     */
    bool intersects(const Interval &x) const;

    /**
     * \brief True iff *this and \a x intersect and the intersection has a non-null volume.
     *
     * Equivalently, some interior points (of this or x) must belong to the intersection.
     */
    bool overlaps(const Interval &x) const;

    /**
     * \brief True iff *this and \a x do not intersect.
     *
     */
    bool is_disjoint(const Interval &x) const;

    /**
     * \brief True iff *this is empty.
     */
    bool is_empty() const;

    /**
     * \brief True iff *this is degenerated.
     *
     * An interval is degenerated if it is of the form [a, a]
     *
     * \note An empty interval is considered here as degenerated. */
    bool is_degenerated() const;

    /**
     * \brief True if one bound of *this is infinite.
     *
     * \note An empty interval is always bounded.
     */
    bool is_unbounded() const;

    /**
     * \brief True iff *this can be bisected into two non-degenerated intervals.
     *
     * Examples of non bisectable intervals are [0,next_float(0)] or [DBL_MAX,+oo).
     */
    bool is_bisectable() const;

    /**
     * \brief Relative Hausdorff distance between *this and x.
     *
     * The relative distance is basically distance(x)/diam(*this).
     * \see #ibex::distance (const ibex::Interval &x1, const ibex::Interval &x2).
     */
    double rel_distance(const Interval& x) const;

    /*
     * \brief The complementary of x.
     */
    int complementary(Interval& c1, Interval& c2, bool compactness = true) const;

    /**
     * \brief x\y
     */
    int diff(const Interval& y, Interval& c1, Interval& c2, bool compactness = true) const;

    /** \brief Return -*this. */
    Interval operator-() const;

    /** \brief Add \a d to *this and return the result.  */
    Interval& operator+=(double d);

    /** \brief Subtract \a d to *this and return the result. */
    Interval& operator-=(double d);

    /** \brief Multiply *this by \a d and return the result. */
    Interval& operator*=(double d);

    /** \brief Divide *this by \a d and return the result. */
    Interval& operator/=(double d);

    /** \brief Add \a x to *this and return the result. */
    Interval& operator+=(const Interval& x);

    /** \brief Subtract \a x to *this and return the result. */
    Interval& operator-=(const Interval& x);

    /** \brief Multiply *this by \a x and return the result. */
    Interval& operator*=(const Interval& x);

    /**
     * \brief Divide *this by \a x and return the result.
     *
     * Does better than *this=*this/x: because calculates
     * the union of *this/x as intermediate result. */
    Interval& operator/=(const Interval& x);

    /**
     * \brief Set this interval to the intersection of itself with the division of two others.
     *
     * \param x - the numerator
     * \param y - the divisor
     * \param out2 - In return, *this and out2 contains the lower and upper part respectively
     * of the division. If the result of the generalized division and intersection
     * is a single interval, out2 is set to the empty interval.
     * \code
     * Interval intv(-10,10);
     * Interval out2;
     * intv.div2_inter(Interval(2,3), Interval(-1,2), out2);
     * cout << intv << " " << out2 << endl;  // will display: [-10,-2] [1,10]
     * \endcode
	 *
     * \return \c true if the intersection is non empty.
     * \note Contrary to the "cset" theory, the result is empty if \a y=[0,0] (whatever \a x is).
     */
    bool div2_inter(const Interval& x, const Interval& y, Interval& out2);

    /**
     * \brief Set this interval to the intersection of itself with the division of two others.
     */
    Interval& div2_inter(const Interval& x, const Interval& y);

    /**
     * \brief Return diam(*this)-diam(x), for x\subseteq *this [deprecated]
     *
     * Deprecated. Kept for compatibility with ibex 1.xx.
     *
     * \pre \a x must be included in this interval.
     * \note The result may be +oo (if the set difference is infinite).
     * \note An empty interval is considered here to have a null diamater (as a degenerated interval). <br>
     * If either \a x or this interval is empty, then the method returns the diameter of this interval
     * (which is 0 if the latter is empty).
     */
    double delta(const Interval& x) const;

    /**
     * \brief Compute the ratio of the diameter to #delta(x) [deprecated].
     *
     * Deprecated. Kept for compatibility with ibex 1.xx.
     *
     * \pre \a x must be included in this interval.
     * \note An empty interval is considered to have a null diamater (as a degenerated interval). <br>
     * <ul><li>If either \a x or this interval is empty, then
     * <ul><li>the method returns 1 (100% of reduction) if this diameter is not null,
     *     <li>0 otherwise (as if 0/0=0).</ul>
     * <li>As a pure convention, the method returns \c 1 if one bound of this interval is infinite and the corresponding bound of \a x
     * is not (in particular if this interval is unbounded and \a x not). </ul>
     */
    double ratiodelta(const Interval& x) const;

    /**
     * \brief Bisect *this into two subintervals.
     *
     * \param ratio - says where to split (0.5=middle)
     * \pre is_bisectable() must be true.
     * \pre 0<ratio<1.
     */
    std::pair<Interval,Interval> bisect(double ratio=0.5) const;

    /** \brief pi
     *  Deprecated. Use pi().
     */
    static const Interval PI;

     /** \brief pi. */
    static const Interval& pi();

    /** \brief 2*pi.
     *  Deprecated. Use two_pi().
     */
    static const Interval TWO_PI;

    /** \brief 2*pi. */
    static const Interval& two_pi();

    /** \brief pi/2.
     *  Deprecated. Use half_pi().
     */
    static const Interval HALF_PI;

    /** \brief pi/2. */
    static const Interval& half_pi();

    /** \brief the empty interval.
     *  Deprecated. Use empty_set().
     */
    static const Interval EMPTY_SET;

    /** \brief the empty interval. */
    static const Interval& empty_set();

    /** \brief (-oo,oo).
     *  Deprecated. Use all_reals().
     */
    static const Interval ALL_REALS;

    /** \brief (-oo,oo). */
    static const Interval& all_reals();

    /** \brief [0,0].
     *  Deprecated. Use zero().
     */
    static const Interval ZERO;

    /** \brief [0,0]. */
    static const Interval& zero();

    /** \brief [1,1].
     *  Deprecated. Use one().
     */
    static const Interval ONE;

    /** \brief [1,1]. */
    static const Interval& one();

    /** \brief [0,+oo).
     *  Deprecated. Use pos_reals().
     */
    static const Interval POS_REALS;

    /** \brief [0,+oo). */
    static const Interval& pos_reals();

    /** \brief (-oo,0].
     *  Deprecated. Use neg_reals().
     */
    static const Interval NEG_REALS;

    /** \brief (-oo,0]. */
    static const Interval& neg_reals();

//    friend class IntervalVector;

    typedef Interval SCALAR;
    typedef IntervalVector VECTOR;
    typedef IntervalMatrix MATRIX;

    /**
     * \brief Cast the interval to an expression
     */
    operator const ExprConstant&() const;

    /* \brief Constructor from an interval [x] of wrapped type. */
    Interval(const interval_type_wrapper& x);
    /* \brief Assign to an interval from a wrapper type. */
    Interval& operator=(const interval_type_wrapper& x);

    interval_type_wrapper itv;
};

/** \ingroup arithmetic */
/*@{*/

/** \brief Stream out \a x. */
std::ostream& operator<<(std::ostream& os, const Interval& x);

/** \brief \f$[x]_1\cap [x]_2\f$.
 * \return Interval::EMPTY if the intersection is empty. */
Interval operator&(const Interval& x1, const Interval& x2);

/** \brief \f$\square([x]_1\cup [x]_2)\f$. */
Interval operator|(const Interval& x1, const Interval& x2);

/** \brief [x]+d. */
Interval operator+(const Interval& x, double d);

/** \brief [x]-d. */
Interval operator-(const Interval& x, double d);

/** \brief [x]*d. */
Interval operator*(const Interval& x, double d);

/** \brief [x]/d. */
Interval operator/(const Interval& x, double d);

/** \brief Generalized division
 *
 * The result of x/y is stored in out1 and out2. These intervals contain the lower and upper part respectively
 * of the division. If the result of the generalized division is a single interval, out2
 * is set to the empty interval. If the result is the empty set, both out1 and out2 are set to the
 * empty interval.
 *
 * \note Contrary to the "cset" theory, the result is empty if \a y=[0,0] (whatever \a x is).
 */
void div2(const Interval& x, const Interval& y, Interval& out1, Interval& out2);

/** \brief d+[x]. */
Interval operator+(double d,const Interval& x);

/** \brief d-[x]. */
Interval operator-(double d, const Interval& x);

/** \brief d*[x]. */
Interval operator*(double d, const Interval& x);

/** \brief d/[x]. */
Interval operator/(double d, const Interval& x);

/** \brief [x]_1+[x]_2. */
Interval operator+(const Interval& x1, const Interval& x2);

/** \brief [x]_1-[x]_2. */
Interval operator-(const Interval& x1, const Interval& x2);

/** \brief [x]_1*[x]_2. */
Interval operator*(const Interval& x1, const Interval& x2);

/** \brief [x]_1/[x]_2. */
Interval operator/(const Interval& x1, const Interval& x2);

/** \brief Hausdorff distance of [x]_1 and [x]_2. */
double distance(const Interval &x1, const Interval &x2);

/** \brief [x]^2 */
Interval sqr(const Interval& x);

/** \brief sqrt{[x]} */
Interval sqrt(const Interval& x);

/** \brief [x]^n. */
Interval pow(const Interval& x, int n);

/** \brief [x]^d. */
Interval pow(const Interval& x, double d);

/** \brief [x]^[y]. */
Interval pow(const Interval &x, const Interval &y);

/** \brief n^{th} root of [x]. */
Interval root(const Interval& x, int n);

/** \brief exp([x]). */
Interval exp(const Interval& x);

/** \brief log([x]). */
Interval log(const Interval& x);

/** \brief cos([x]). */
Interval cos(const Interval& x);

/** \brief sin([x]). */
Interval sin(const Interval& x);

/** \brief tan([x]). */
Interval tan(const Interval& x);

/** \brief acos([x]). */
Interval acos(const Interval& x);

/** \brief asin([x]). */
Interval asin(const Interval& x);

/** \brief atan([x]). */
Interval atan(const Interval& x);

/** \brief atan2([y],[x]). */
Interval atan2(const Interval& y, const Interval& x);

/** \brief cosh([x]). */
Interval cosh(const Interval& x);

/** \brief sinh([x]). */
Interval sinh(const Interval& x);

/** \brief tanh([x]). */
Interval tanh(const Interval& x);

/** \brief acosh([x]). */
Interval acosh(const Interval& x);

/** \brief asinh([x]). */
Interval asinh(const Interval& x);

/** \brief atanh([x]). */
Interval atanh(const Interval& x);

/** \brief \f$abs([x]) = \{|x|, x\in[x]\}.\f$. */
Interval abs(const Interval &x);

/** \brief Maximum of two intervals.
 *
 * Return \f$\max([x],[y]) = \{ \max\{x,y\}, x\in[x], y\in[y] \}\f$. */
Interval max(const Interval& x, const Interval& y);

/** \brief Minimum of two intervals.
 *
 *  Return \f$\min([x],[y]) = \{ \min\{x,y\}, x\in[x], y\in[y] \}\f$. */
Interval min(const Interval& x, const Interval& y);

/** \brief Sign of [x].
 *
 *  Return \f$sign([x]) = hull \{ sign{x}, x\in[x] \}\f$.
 * \remark By convention, \f$ 0\in[x] \Longrightarrow sign[x]=[-1,1]\f$. */
Interval sign(const Interval& x);

/** \brief Chi of [a], [b] and [c].
 *
 *  Return \f$chi([a],[b],[c]) = 0.5*(1-\sign([a]))*[b] + 0.5*(\sign([a])+1)*[c]. \f$
 * \remark  chi([a],[b],[c]) =[b] if [a]<=0, [c] if [a]>0, hull \{[b], [c]\} else.  */
Interval chi(const Interval& a, const Interval& b, const Interval& c);

/** \brief Return the largest integer interval included in x. */
Interval integer(const Interval& x);

/** \brief Floor of [x]. */
Interval floor(const Interval& x);

/** \brief Ceil of [x]. */
Interval ceil(const Interval& x);

/**
 * \brief The saw function is x->x-round(x)
 *
 * It allows to enforce integer constraint
 * by writing saw(x)=0.
 */
Interval saw(const Interval& x);

/** \brief Projection of y=x_1+x_2.
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=x_1+x_2\}\f$. */
bool bwd_add(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=x_1-x_2.
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=x_1-x_2\}\f$. */
bool bwd_sub(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=x_1*x_2.
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=x_1\times x_2\}\f$. */
bool bwd_mul(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=x_1/x_2.
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=x_1/x_2\}\f$. */
bool bwd_div(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=x^2.
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=x^2 \}\f$. */
bool bwd_sqr(const Interval& y, Interval& x);

/** \brief Projection of y=sqrt{x}.
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=sqrt{x} \}\f$. */
bool bwd_sqrt(const Interval& y, Interval& x);

/** \brief Projection of y=x^n.
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=x^n \}\f$. */
bool bwd_pow(const Interval& y, int n, Interval& x);

/** \brief Projection of y=x_1^{x_2}.
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=x_1^{x_2}\}\f$. */
bool bwd_pow(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of the \f$y=x^{\frac{1}{n}}\f$.
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=x^{\frac{1}{n}} \}\f$. */
bool bwd_root(const Interval& y, int n, Interval& x);

/** \brief Projection of y=exp(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\exp(x) \}\f$. */
bool bwd_exp(const Interval& y,  Interval& x);

/** \brief Projection of y=log(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\log(x) \}\f$. */
bool bwd_log(const Interval& y,  Interval& x);

/** \brief Projection of y=cos(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\cos(x) \}\f$. */
bool bwd_cos(const Interval& y,  Interval& x);

/** \brief Projection of y=sin(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\sin(x) \}\f$. */
bool bwd_sin(const Interval& y,  Interval& x);

/** \brief Projection of y=tan(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\tan(x) \}\f$. */
bool bwd_tan(const Interval& y,  Interval& x);

/** \brief Projection of y=acos(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\arccos(x) \}\f$. */
bool bwd_acos(const Interval& y,  Interval& x);

/** \brief Projection of y=asin(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\arcsin(x) \}\f$. */
bool bwd_asin(const Interval& y,  Interval& x);

/** \brief Projection of y=atan(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\arctan(x) \}\f$. */
bool bwd_atan(const Interval& y,  Interval& x);

/** \brief Projection of y=atan2(x_1,x_2).
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=atan2(x_1,x_2)\f$. */
bool bwd_atan2(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=cosh(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\cosh(x) \}\f$. */
bool bwd_cosh(const Interval& y,  Interval& x);

/** \brief Projection of y=sinh(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\sinh(x) \}\f$. */
bool bwd_sinh(const Interval& y,  Interval& x);

/** \brief Projection of y=tanh(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=\tanh(x) \}\f$. */
bool bwd_tanh(const Interval& y,  Interval& x);

/** \brief Projection of y=acosh(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=arccosh(x) \}\f$. */
bool bwd_acosh(const Interval& y,  Interval& x);

/** \brief Projection of y=asinh(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=arcsinh(x) \}\f$. */
bool bwd_asinh(const Interval& y,  Interval& x);

/** \brief Projection of y=atanh(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=arctanh(x) \}\f$. */
bool bwd_atanh(const Interval& y,  Interval& x);

/** \brief Projection of y=|x|.
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=|x| \}\f$. */
bool bwd_abs(const Interval& y,  Interval& x);

/** \brief Projection of y=max(x_1,x_2).
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=\max(x_1,x_2)\}\f$. */
bool bwd_max(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=min(x_1,x_2).
 *
 * Set ([x]_1,[x]_2) to \f$([x]_1,[x]_2])\cap\{ (x_1,x_2)\in [x]_1\times[x]_2 \ | \ \exists y\in[y],\ y=\min(x_1,x_2)\}\f$. */
bool bwd_min(const Interval& y, Interval& x1, Interval& x2);

/** \brief Projection of y=sign(x).
 *
 * Set [x] to \f$[x]\cap \{ x\in [x] \exists y\in [y], \quad y=sign(x) \}\f$. */
bool bwd_sign(const Interval& y, Interval& x);

/** \brief Projection of f=chi(a,b,c). */
bool bwd_chi(const Interval& f, Interval& a, Interval& b, Interval& c);

/** \brief Contract x w.r.t. y=floor(x). */
bool bwd_floor(const Interval& y, Interval& x);

/* \brief Contract x w.r.t. y=ceil(x). */
bool bwd_ceil(const Interval& y, Interval& x);

/* \brief Contract x w.r.t. y=saw(x). */
bool bwd_saw(const Interval& y, Interval& x);

/**
 * \brief Contract x and y w.r.t. the fact that they are equivalent modulo the period p.
 */
bool bwd_imod(Interval& x, Interval& y, const double& p);

} // end namespace ibex

#include "ibex_IntervalLibWrapper.inl"

/*@}*/

/*============================================ features with common implementation ============================================ */

namespace ibex {

// the following functions are
// introduced to allow genericity
inline bool ___is_empty(double)             { return false; }
inline bool ___is_empty(const Interval& x)  { return x.is_empty(); }
inline void ___set_empty(double)            { }
inline void ___set_empty(Interval& x)       { x.set_empty(); }
inline double ___abs(double x)              { return fabs(x);}
inline Interval ___abs(const Interval& x)   { return abs(x);}
inline double ___mag(double x)              { return fabs(x);}
inline double ___mag(const Interval& x)     { return x.mag();}

inline Interval::Interval() : itv(NEG_INFINITY, POS_INFINITY) {

}

inline Interval::Interval(double a, double b) : itv(a,b) {
	if (a==POS_INFINITY || b==NEG_INFINITY || a>b) *this=EMPTY_SET;
}

inline Interval::Interval(double a) : itv(a,a) {
	if (a==NEG_INFINITY || a==POS_INFINITY) *this=EMPTY_SET;
}

inline Interval::Interval(std::array<double, 1> array): Interval(array[0]) {}

inline Interval::Interval(std::array<double, 2> array): Interval(array[0], array[1]) {}

inline bool Interval::operator==(const Interval& x) const {
	return (is_empty() && x.is_empty()) || (lb()==x.lb() && ub()==x.ub());
}

inline Interval& Interval::operator=(double x) {
	if (x==NEG_INFINITY || x==POS_INFINITY)
		*this=EMPTY_SET;
	else
		itv = x;
	return *this;
}

inline Interval& Interval::operator=(const Interval& x) {
	itv = x.itv;
	return *this;
}

inline Interval& Interval::inflate(double radd) {
	(*this) += Interval(-radd,radd);
	return *this;
}

inline Interval& Interval::inflate(double delta, double chi) {
	double m=mid();
	(*this) = m + delta*(*this-m)+Interval(-chi,chi);
	return *this;
}


inline bool Interval::operator!=(const Interval& x) const {
	return !(*this==x);
}

inline double Interval::rad() const {
	if (is_empty()) return 0;
	else if (is_unbounded()) return POS_INFINITY;
	else {
		double t = mid();
		double t1 =(t-*this).ub();
		double t2= (*this-t).ub();
		return (t1>t2) ? t1 : t2;
	}
}

inline bool Interval::is_bisectable() const {
	if (is_empty()) return false;
	double m=mid();
	return (lb()<m && m<ub());
}


inline double Interval::rel_distance(const Interval& x) const {
	  double d=distance(*this,x);
	  if (d==POS_INFINITY) return 1;
	  double D=diam();
	  return (D==0 || D==POS_INFINITY) ? 0.0 : (d/D); // if diam(this)=infinity here, necessarily d=0
}

inline double distance(const Interval &x1, const Interval &x2) {

    if (x1.is_empty()) return x2.rad();

    if (x2.is_empty()) return x1.rad();

    if (x1.lb()==NEG_INFINITY) {
    	if (x2.lb()!=NEG_INFINITY)
    		return POS_INFINITY;
    	else if (x1.ub()==POS_INFINITY) {
    		if (x2.ub()==POS_INFINITY) return 0.0;
    		else return POS_INFINITY;
    	}
    	else if (x2.ub()==POS_INFINITY) return POS_INFINITY;
    	else return fabs(x1.ub()-x2.ub());
    }
    else if (x1.ub()==POS_INFINITY) {
    	if (x2.ub()!=POS_INFINITY)
    		return POS_INFINITY;
    	else if (x2.lb()==NEG_INFINITY)
    		return POS_INFINITY;
    	else return fabs(x1.lb()-x2.lb());
    }
    else if (x2.is_unbounded())
    	return POS_INFINITY;
    else
    	return _interval_distance_wrapper (x1.itv, x2.itv);
}

inline Interval sign(const Interval& x) {
	return x.ub()<0 ? Interval(-1,-1) : x.lb()>0 ? Interval(1,1) : Interval(-1,1);
}

inline Interval chi(const Interval& a, const Interval& b, const Interval& c){
	if (a.ub()<=0) {
		return b;
	} else if (a.lb()>0) {
		return c;
	} else {
		return b|c;
	}
}

inline Interval atan2(const Interval& y, const Interval& x) {

	if (y.is_empty() || x.is_empty()) return Interval::empty_set();

	// we handle the special case x=[0,0] separately
	else if (x==Interval::zero()) {
		if (y.lb()>=0)
			if (y.ub()==0) return Interval::empty_set(); // atan2(0,0) is undefined.
			else return Interval::half_pi();
		else if (y.ub()<=0) return -Interval::half_pi();
		else return Interval(-1,1)*Interval::half_pi();
	}

	else if (x.lb()>=0) {
		return atan(y/x); // now, x.ub()>0 -> atan does not give an empty set
	}
	else if (x.ub()<=0) {
		if (y.lb()>=0)
			return atan(y/x)+Interval::pi(); // x.lb()<0
		else if (y.ub()<0)
			return atan(y/x)-Interval::pi();
		else
			return Interval(-1,1)*Interval::pi();
	} else {
		if (y.lb()>=0)
			return atan(y/x.ub()) | (atan(y/x.lb()) + Interval::pi());
		else if (y.ub()<=0){
			if(x.lb()!=NEG_INFINITY){
				if(x.ub()!=POS_INFINITY){
					return (atan(y/x.lb())-Interval::pi()) | atan(y/x.ub());
				}
				else
					return (atan(y/x.lb())-Interval::pi()) | Interval::zero();
			}
			else{
				if(x.ub()!=POS_INFINITY)
					return (-Interval::pi()) | atan(y/x.ub());
				else
					return -Interval::pi() | Interval::zero();
			}
		}
		else
			return Interval(-1,1)*Interval::pi();
	}
}

inline bool bwd_add(const Interval& y, Interval& x1, Interval& x2) {
	if ((x1 &= y-x2).is_empty()) { x2.set_empty(); return false; }
	if ((x2 &= y-x1).is_empty()) { x1.set_empty(); return false; }
	return true;
}

inline bool bwd_sub(const Interval& y, Interval& x1, Interval& x2) {
	if ((x1 &= y+x2).is_empty()) { x2.set_empty(); return false; }
	if ((x2 &= x1-y).is_empty()) { x1.set_empty(); return false; }
	return true;
}

inline bool bwd_div(const Interval& y, Interval& x1, Interval& x2) {
	if ((x1 &= y*x2).is_empty()) { x2.set_empty(); return false; }
	Interval tmp=y;
	bwd_mul(x1, tmp, x2);
	if (x2.is_empty()) { x1.set_empty(); return false; }
	return true;
}

inline bool bwd_sqrt(const Interval& y, Interval& x) {
	if (y.is_empty() || y.ub()<0) {
		x.set_empty();
	} else if (y.lb()<0) {
		x &= sqr(Interval(0,y.ub()));
	} else  {
		x &= sqr(y);
	}
	return !x.is_empty();
}

inline bool bwd_root(const Interval& y, int n, Interval& x) {
	x &= pow(y,n);
	return !x.is_empty();
}

inline bool bwd_exp(const Interval& y,  Interval& x) {
	x &= log(y);
	return !x.is_empty();
}

inline bool bwd_log(const Interval& y,  Interval& x) {
	x &= exp(y);
	return !x.is_empty();
}

inline bool bwd_acos(const Interval& y,  Interval& x) {
	x &= cos(y);
	return !x.is_empty();
}

inline bool bwd_asin(const Interval& y,  Interval& x) {
	x &= sin(y);
	return !x.is_empty();
}

inline bool bwd_atan(const Interval& y,  Interval& x) {

	if (y.is_empty()) {
		x.set_empty();
		return false;
	}

	// Note: if y.ub>pi/2 or y.lb<-pi/2, tan(y) gives (-oo,oo).
	// so the implementation is not as simple as x &= tan(y).

	Interval z=y;
	double pi2l=(Interval::pi()/2).lb();
	double pi2u=(Interval::pi()/2).ub();

	if (z.ub()>=pi2l) // not pi2u. See comments below.
		if (z.lb()>=pi2u)
			x.set_empty();
		else {
			if (z.lb()>-pi2l) {
				// note 1: tan(z^-) can give an interval (-oo,+oo) if
				// z^- is close to -pi/2. Even in this case we keep the
				// lower bound -oo.
				//
				// note 2: if we had used z.lb()<-pi2u (with pi2u>pi/2)
				// instead of z.lb()<-pi2l, it may be possible, in theory,
				// that the calculated lower bound is a high value close to +oo, which would be incorrect.
				//
				// note 3: if z.lb() is close to pi/2, the lower bound of tan(z.lb()) can be -oo. There
				// is nothing we can do about it (the lower bound cannot be evaluated in this case)
				//
				// note 4: tan(z.lb()) cannot be an empty set since z.lb() cannot be exactly pi/2.
				x &= Interval(tan(Interval(z.lb())).lb(),POS_INFINITY);
			}
			// else do nothing
		}
	else
		if (z.ub()<=-pi2u)
			x.set_empty();
		else if (z.lb()<-pi2l)
			// Same comments as above.
			x &= Interval(NEG_INFINITY,tan(Interval(z.ub())).ub());
		else
			x &= Interval(tan(Interval(z.lb())).lb(),
					tan(Interval(z.ub())).ub());

	return !x.is_empty();
}


inline bool bwd_acosh(const Interval& y,  Interval& x) {
	if (y.is_empty() || y.ub()<0.0) {
		x.set_empty(); return false;
	}
	else {
		x &= cosh(Interval(y.lb()<0?0:y.lb(),y.ub()));
		return !x.is_empty();
	}
}

inline bool bwd_asinh(const Interval& y,  Interval& x) {
	x &= sinh(y);
	return !x.is_empty();
}

inline bool bwd_atanh(const Interval& y,  Interval& x) {
	x &= tanh(y);
	return !x.is_empty();
}

inline bool bwd_max(const Interval& y, Interval& x1, Interval& x2) {

	if (y.is_empty()) {
		x1.set_empty();
		x2.set_empty();
		return false;
	}

	/* ---- Disjoint intervals ---- */
	if (x2.lb()>x1.ub() || y.lb()>x1.ub()) {
		/* then, max(x,x2) is necessarily x2 */
		if ((x2 &= y).is_empty()) { x1.set_empty(); return false;}
		else return true;
	} else if (x1.lb()>x2.ub() || y.lb()>x2.ub()) {
		if ((x1 &= y).is_empty()) { x2.set_empty(); return false;}
		else return true;
	}
	/*------------------------------*/

	if (y.ub()<x1.lb() || y.ub()<x2.lb()) {
		x1.set_empty();
		x2.set_empty();
		return false; // inconsistency
	}

	/* At this point, x1, x2 and y all mutually intersect. */
	if (x1.ub()>y.ub())
		x1=Interval(x1.lb(),y.ub());
	if (x2.ub()>y.ub())
		x2=Interval(x2.lb(),y.ub());

	return true;
}

inline bool bwd_min(const Interval& y, Interval& x1, Interval& x2) {

	Interval mx1=-x1;
	Interval mx2=-x2;

	if (!bwd_max(-y,mx1,mx2)) {
		x1.set_empty();
		x2.set_empty();
		return false;
	}

	x1=-mx1;
	x2=-mx2;
	return true;
}

inline bool bwd_sign(const Interval& y,  Interval& x) {

	if (y.is_empty()) {
		x.set_empty();
		return false;
	}

	if(y.lb()>0)
		x &= Interval::pos_reals();
	else if(y.ub()<0)
		x &= Interval::neg_reals();
	return !x.is_empty();
}


// probably not the most efficient implementation ever...
inline bool bwd_atan2(const Interval& theta, Interval& y, Interval& x) {

	// half-plane x>0
	Interval theta_xpos = theta  & Interval(-1,1)*Interval::half_pi();
	// upper left quadrant
	Interval theta_xneg_ypos = theta  & (Interval::half_pi() | Interval::pi());
	// lower left quadrant
	Interval theta_xneg_yneg = theta  & -(Interval::half_pi() | Interval::pi());

	Interval xres= Interval::empty_set();
	Interval yres= Interval::empty_set();

	if (!theta_xpos.is_empty()) {

		Interval xpos=x & Interval::pos_reals();
		Interval yall=y;

		if (theta_xneg_ypos.is_empty() || theta_xneg_yneg.is_empty()) { // otherwise, all is valid (useless to contract)
			Interval z=yall/xpos;
			bwd_atan(theta_xpos,z);
			bwd_div(z,yall,xpos);

		}
		xres |= xpos;
		yres |= yall;

		// because atan is not defined for pi/2 and -pi/2
		if (theta_xpos.lb()>=Interval::half_pi().lb()) {
			xres |= (x & Interval::zero());
			yres |= (y & Interval::pos_reals());
		} else if (theta_xpos.ub()<=-Interval::half_pi().lb()) {
			xres |= (x & Interval::zero());
			yres |= (y & Interval::neg_reals());
		}
	}

	if (!theta_xneg_ypos.is_empty()) {
		Interval xneg=x & Interval::neg_reals();
		Interval ypos=y & Interval::pos_reals();
		Interval z=ypos/xneg;
		bwd_atan(theta_xneg_ypos - Interval::pi(),z);
		bwd_div(z,ypos,xneg);
		xres |= xneg;
		yres |= ypos;	}

	if (!theta_xneg_yneg.is_empty()) {
		Interval xneg=x & Interval::neg_reals();
		Interval yneg=y & Interval::neg_reals();
		Interval z=yneg/xneg;
		bwd_atan(theta_xneg_yneg + Interval::pi(),z);
		bwd_div(z,yneg,xneg);
		xres |= xneg;
		yres |= yneg;
	}

	x=xres;
	y=yres;

	return !x.is_empty();
}


//inline bool bwd_atan2(const Interval& theta, Interval& y, Interval& x) {
//
//	if (theta.is_empty()) {
//		x.set_empty(); y.set_empty();
//		return false;
//	}
//
//    //Lower half of upper right quadrant
//    if(theta.is_subset(Interval(0,M_PI/4.)))
//    {
//        x = (x&Interval::pos_reals()) & ( (y&Interval::pos_reals()) * (1/tan(theta&Interval(0,M_PI/4.))) );
//        y = (y&Interval::pos_reals()) & ( (x&Interval::pos_reals()) * (tan(theta&Interval(0,M_PI/4.))) );
//    }
//
//    //Upper half of upper right quadrant
//    else if(theta.is_subset(Interval(M_PI/4.,M_PI/2.)))
//    {
//        bwd_atan2(M_PI/2.-theta,x,y);
//    }
//
//    //Upper right quadrant
//    else if(theta.is_subset(Interval(0,M_PI/2)))
//    {
//        Interval x1=x; Interval y1=y; Interval x2=x; Interval y2=y;
//        Interval theta1(theta.lb(),M_PI/4.), theta2(M_PI/4., theta.ub());
//        bwd_atan2(theta1,y1,x1);
//        bwd_atan2(theta2,y2,x2);
//        x=x1|x2; y=y1|y2;
//    }
//
//    //Upper left quadrant
//    else if(theta.is_subset(Interval(M_PI/2,M_PI)))
//    {
//        Interval x2=-x;
//        Interval theta2=M_PI-theta;
//        bwd_atan2(theta2,y,x2);
//        x=-x2;
//    }
//
//    //Lower left quadrant
//    else if(theta.is_subset(Interval(M_PI,3*M_PI/2)))
//    {
//        Interval y2=-y;
//        Interval x2=-x;
//        Interval theta2=theta-M_PI;
//        bwd_atan2(theta2,y2,x2);
//        x=-x2; y=-y2;
//    }
//
//    //Lower right quadrant
//    else if(theta.is_subset(Interval(3*M_PI/2,2*M_PI)))
//    {
//        Interval x2=x;
//        Interval y2=-y;
//        Interval theta2=2*M_PI-theta;
//        bwd_atan2(theta2,y2,x2);
//        x=x2; y=-y2;
//    }
//
//    //Upper bissection
//    else if(theta.is_subset(Interval(0,M_PI))) {
//        Interval theta1(theta.lb(),M_PI/2), theta2(M_PI/2, theta.ub());
//        Interval x1=x; Interval y1=y; Interval x2=x; Interval y2=y;
//        if(theta.lb() != M_PI/2) bwd_atan2(theta1,y1,x1);
//        if(theta.ub() != M_PI/2) bwd_atan2(theta2,y2,x2);
//        x=x1|x2; y=y1|y2;
//    }
//
//    //Left bissection
//    else if(theta.is_subset(Interval(0,3*M_PI/2))) {
//        Interval theta1(theta.lb(),M_PI), theta2(M_PI, theta.ub());
//        Interval x1=x; Interval y1=y; Interval x2=x; Interval y2=y;
//        if(theta.lb() != M_PI) bwd_atan2(theta1,y1,x1);
//        if(theta.ub() != M_PI) bwd_atan2(theta2,y2,x2);
//        x=x1|x2; y=y1|y2;
//    }
//
//    //Lower bissection
//    else if(theta.is_subset(Interval(0,2*M_PI))) {
//        Interval theta1(theta.lb(),3*M_PI/2), theta2(3*M_PI/2, theta.ub());
//        Interval x1=x; Interval y1=y; Interval x2=x; Interval y2=y;
//        if(theta.lb() != 3*M_PI/2) bwd_atan2(theta1,y1,x1);
//        if(theta.ub() != 3*M_PI/2) bwd_atan2(theta2,y2,x2);
//        x=x1|x2; y=y1|y2;
//    }
//
//    //Theta diameter greater than 2PI then theta will be considered as [0,2PI]
//    else if (theta.diam() > 2*M_PI)
//    {
//        Interval theta1(0,2*M_PI);
//        bwd_atan2(theta1,x,y);
//    }
//
//    // Modulo
//    else
//    {
//        // We separate the intervals into 4 cases as imod does not support union.
//        Interval theta1(0,M_PI/2.), theta2(M_PI/2., M_PI), theta3(M_PI, 3*M_PI/2.), theta4(3*M_PI/2., 2*M_PI);
//        Interval thetaTmp(theta);
//        bwd_imod(thetaTmp,theta1,2*M_PI);
//        thetaTmp=theta;
//        bwd_imod(thetaTmp,theta2,2*M_PI);
//        thetaTmp=theta;
//        bwd_imod(thetaTmp,theta3,2*M_PI);
//        thetaTmp=theta;
//        bwd_imod(thetaTmp,theta4,2*M_PI);
//        Interval x1=x; Interval y1=y;
//        bwd_atan2(theta1,y1,x1); // first quadrant
//        Interval x2=x; Interval y2=y;
//        bwd_atan2(theta2,y2,x2); // second quadrant
//        Interval x3=x; Interval y3=y;
//        bwd_atan2(theta3,y3,x3); // third quadrant
//        Interval x4=x; Interval y4=y;
//        bwd_atan2(theta4,y4,x4); // fourth quadrant
//        x=(x1|x2)|(x3|x4); y=(y1|y2)|(y3|y4);
//        // not_implemented("bwd_atan2 not implemented yet for theta outside [0,2*PI].");
//    }
//
//
//    return !x.is_empty() || !y.is_empty();
//}

inline bool bwd_chi(const Interval& f, Interval& a, Interval& b, Interval& c){
	if      (a.ub()<=0) {if ((b &= f).is_empty()) { a.set_empty(); c.set_empty(); return false; } }
	else if (a.lb()>0)  {if ((c &= f).is_empty()) { a.set_empty(); b.set_empty(); return false; } }

	if (f.is_disjoint(b)) {
		if ((a &= Interval::pos_reals()).is_empty()) { b.set_empty(); c.set_empty(); return false; }
		if ((c &= f).is_empty()) { a.set_empty(); b.set_empty(); return false; }
	}
	if (f.is_disjoint(c)) {
		if ((a &= Interval::neg_reals()).is_empty()) { b.set_empty(); c.set_empty(); return false; }
		if ((b &= f).is_empty()) { a.set_empty(); c.set_empty(); return false; }
	}
	return true;
}


inline bool bwd_floor(const Interval& y, Interval& x) {
	if (y.is_empty()) {
		x.set_empty();
		return false;
	} else {
		double r=std::floor(y.ub());
		double l=std::ceil(y.lb());
		if (r<l) {
			x.set_empty();
			return false;
		}
		else {
			x &= Interval(l,r) + Interval(0,1);
			return !x.is_empty();
		}
	}
}

inline bool bwd_ceil(const Interval& y, Interval& x) {
	if (y.is_empty()) {
		x.set_empty();
		return false;
	} else {
		double r=std::floor(y.ub());
		double l=std::ceil(y.lb());
		if (r<l) {
			x.set_empty();
			return false;
		}
		else {
			x &= Interval(l,r) + Interval(-1,0);
			return !x.is_empty();
		}
	}
}

// Implements interval modulo with double period:  x = y mod(p)
inline bool bwd_imod(Interval& x, Interval& y, const double& p) {
    if (p <= 0.)
    {
        ibex_error("Modulo needs a strictly positive period p.");
        return false;
    }
    if (y.diam()>p || x.diam()>p)
        return false;
    Interval r = (x-y)/p;
    Interval ir = integer(r);
    if (ir.is_empty()) // additional protection because an empty interval is considered degenerated.
    {
        x.set_empty(); y.set_empty();
        return false;
    }
    if (ir.is_degenerated())
        bwd_sub(ir*p,x,y);
    else if (ir.diam()==1.)
    {
        double ir1 = ir.lb();
        double ir2 = ir.ub();
        Interval x1 = x; Interval x2 = x;
        Interval y1 = y; Interval y2 = y;
        bwd_sub(Interval(ir1*p),x1,y1);
        bwd_sub(Interval(ir2*p),x2,y2);
        x = x1 | x2;
        y = y1 | y2;
    }
    else
    {
        ibex_error("Modulo diameter error.");
        return false;
    }

    return true;
}

} // end namespace ibex

#endif // _IBEX_INTERVAL_H_
